Journal of The Royal Society Interface

Kinetic parameters in engineered biomolecular circuits are typically known only within intervals and are not precisely tunable; controller synthesis is therefore naturally posed over sets of gains rather than at a single nominal design. We formulate robust controller tuning under parametric plant uncertainty as a universally quantified set-inversion problem: given a plant-parameter box, compute the controller gains that stabilize every plant in the box. Robust stability is decided by Hurwitz sign conditions, which are bounded over the joint gain-parameter domain using validated range enclosures, yielding machine-checkable certificates of local exponential stability of the deterministic equilibrium. Numerical experiments show that the certified region depends strongly on the range enclosure. For an antithetic integral feedback (AIF) controller acting on a two-species gene-expression cascade, a tensor-product Bernstein enclosure certifies 83.2% of the design plane at bisection tolerance epsilon=0.05 with 324 terminal boxes, whereas four baseline enclosures certify between 50.9% and 58.9% on the same paving. For a fifth-order loop, the certified fraction increases from 34.4% to 73.8%, and the certified and excluded terminal boxes confine the robustly stabilizing area to a window of about 2.1 percentage points. Finally, for a refining variant of the algorithm in which the plant box is refined alongside the design box, the undecided area satisfies |U(epsilon)|=O(epsilon) under a non-degeneracy assumption on the stability boundary.

Forecasting and causal inference in complex, non-linear dynamical systems are fundamentally challenging due to noisy and high-dimensional time-series data. While Empirical Dynamic Modeling addresses this via state-space reconstruction in Euclidean space, it suffers from the "curse of dimensionality," distorting distance measurements as dimensions increase. This study proposes Forest-Weighted S-map (FORWS) and Forest-Weighted Causal Inference (FORWC), utilizing adaptive "forest weights" derived from random forest ensembles instead of Euclidean metrics. Using simulated and empirical ecological datasets, I demonstrate that these methods exhibit comparable or improved forecasting skill relative to conventional tools, while offering substantial resilience against dynamic process noise and mitigating the curse of dimensionality. Notably, FORWC enables multimodal causal inference, successfully extracting directed interactions between high-dimensional acoustic vectors and scalar temperature data monitored in a honeybee hive. By bridging chaotic physics with machine learning, this framework accurately extracts intrinsic manifold topology, providing a robust paradigm for uncovering causal networks in complex real-world environments.

Interception of moving objects is a universal pattern of animal behavior with many industrial applications. Inspired by the behavior of the characid fish Brycon guatemalensis, which catches fruits falling onto the surface of a river, we study one-shot (i.e., open-loop) interceptions using both stochastic optimal control theory and reinforcement learning. In the first instance, an agent actuates a control (wait time before launch and heading angle) to meet its target, given a noisy estimate of the target initial state. In a minimal setting, where launch velocity and actuated control are constants, we show that the optimal wait time is attenuated by noise in both the target measurement and control actuation, as well as agent energy conservation. We then derive an upper bound for the heading angle variance which scales as the ratio of the non-vanishing optimal interception error to the initial planar separation. Complementing our control framework, our deep reinforcement learning approach shows that the neural network representing the fish first learns the value function for the interception, then learns how to aim, and finally learns the optimal wait time. Beyond the specific problem at hand, our results may also have applications to problems with little time for feedback or corrections such as perching, landing, and docking.

Behavioral variability reflects differences in internal and external parameters and provides a substrate for natural selection. However, uncovering its structure is challenging because behavior evolves nonlinearly across multiple timescales. Our goal is to reconstruct phenotypic variation from finite behavioral observations without prescribing the scales at which individuals should be compared. To achieve this, we introduce a multiscale dissimilarity that compares dynamics across all statistically resolvable resolutions by constructing reduced transfer operators at progressively finer scales and weighting their differences by finite-sampling uncertainty. This dissimilarity defines a geometry of dynamical phenotypes, which we reconstruct using diffusion maps. We show that the resulting phenotypic coordinates locally follow the parameter directions to which dynamics are most sensitive. We validate our approach in a stochastic model and apply it to bacterial and larval zebrafish behavior. In bacteria, we identify run speed as a major component of inter-individual variability, associated with the chemotaxis protein CheB. In zebrafish, we find that prior prey experience reshapes behavioral dynamics and phenotypic variability, with paramecia experience concentrating fish in a prey-capture-associated region of phenotypic space. Together, these results provide a framework for identifying internal and environmental variables that structure dynamical phenotypes across scales.

Understanding the relationship between weather conditions and infectious disease transmission has important public health implications; nevertheless, the functional form of this relationship often remains poorly understood. Humidity has previously been suggested to have a U-shaped effect on the transmission of respiratory pathogens, with optimal transmission occurring at both low and high levels. Here, to assess the impact of this functional form and of its potential misspecification, we first simulate transmission rates based on the U-shaped relationship and then test whether an exponential function accurately estimates these rates across a range of climates. We further use a transmission model to simulate case count data for a range of hypothetical pathogens, comparing outbreak dynamics under the U-shaped forcing function to those under the best-fit exponential approximation. In our simulations, the model fit for the exponential function is climate-dependent, performing best in very dry and very humid environments. In locations where specific humidity values centre around the vertex of the U-shaped function, the exponential model underestimates the amplitude and frequency of incidence peaks. These findings demonstrate the potential consequences of assuming an incorrect functional form for weather-infectious disease relationships and highlight the need for more careful consideration of this functional form in modelling practice.

Trigger waves are positive feedback-mediated traveling waves that propagate without slowing down or diminishing in amplitude, enabling biochemical coordination over distances where diffusion alone would be insufficient. Apoptosis can spread as trigger waves of caspase activity, as demonstrated in Xenopus eggs and egg extracts. We measured apoptotic wave speeds in extracts at two cell cycle states and found that waves propagated more slowly in M-phase arrested cytostatic factor (CSF) extracts than in interphase extracts (21 versus 27 m min-1). Luther's equation, v=2{surd}D{kappa}, predicts wave speed from the diffusion coefficient D and the positive-feedback doubling rate {kappa}. However, D and {kappa} were comparable between the two extracts and could not explain the speed difference. This discrepancy highlighted a key limitation of the classical equation: it does not account for an activation threshold in trigger wave initiation, and in fact the threshold required to trigger apoptosis was found to differ between cell cycle states. We therefore extended Luther's equation to include such a threshold, obtaining v{propto}{surd}(AD/Uc ), where A is the reaction rate and Uc is the activation threshold. Perturbing apoptotic sensitivity confirmed the predicted threshold-speed relationship. These findings establish the activation threshold as a key determinant of trigger wave speed and provide a broadly applicable scaling law for predicting the spread of switch-like biochemical signals.

We show that the p53-Mdm2 gene interaction can be understood as a form of antithetic feedback, but with mixed feedback loops. In particular, we show that this viewpoint sheds light on the emergence and robustness properties of p53's oscillatory behaviour, as observed experimentally.

From the iconic V-formations of geese to the spectacular murmurations of starlings and communal roosting of millions of bramblings, bird flocks keep sparking wonder and scientific curiosity. In science, common themes include how flocking formation impacts the energy expenditure during migration, how it affects risks of predation, how it relates to information spread on food abundance, and how birds manoeuvre in dense flocks. Here, we report on observations of a very different type of flock formation: migrating terns over open sea flying side by side in a search party manner -- strikingly orthogonal to their heading -- with interbird distances up towards 300 m and flocks of around ten birds covering widths up to one kilometer. Such sparse formations of migrating birds stand in great contrast to the classes of bird flocks previously described in the literature. We present simultaneous field observations and ship-based radar measurements of flock characteristics and argue that this is a form of fly-and-forage migration that may be a common, but overlooked, network foraging strategy of surface feeding birds: while most migrants cannot combine migration and foraging, the situation is different for surface feeders like terns.

Crank length affects the geometry of the bicycle-rider system and influences the motion of the lower limb throughout the pedal cycle. Conventional crank selection has largely relied on rider stature and mechanical power, yet mechanical power and cycling economy remain remarkably insensitive to crank length across a broad range. Yet, elite cyclists and triathletes increasingly adopt shorter cranks, particularly in disciplines where aerodynamic drag strongly affects cycling performance. However, the mechanisms that connect crank length to mechanical power and aerodynamic drag remain unclear. Here we show that crank length couples rider anthropometry, hip kinematics, and aerodynamic drag. Population anthropometry of 6,068 individuals reveals that stature alone poorly predicts personalized crank length, while geometric analysis identifies a nearly linear relationship between crank length and minimum lower hip angle: a 10-mm reduction in crank length increases hip-angle margin by approximately 1.42 degrees with little mechanical penalty. This geometric freedom permits more flexion of the torso and a more aerodynamic posture; the resulting aerodynamic advantage grows rapidly with cycling speed because aerodynamic power scales with velocity cubed. Data from 150 elite cyclists and triathletes support this shift toward shorter cranks, with triathletes adopting shorter relative crank lengths than cyclists. Our results establish crank length as a personalized design variable that connects lower-limb geometry to aerodynamics. These insights enable personalized crank selection, integrated rider-bicycle design, and performance optimization for the increased aerodynamic demands of high-speed cycling.

Blooms of the toxic dinoflagellate Ostreopsis cf. ovata can cause respiratory irritation in beach users, yet their environmental drivers and temporal dynamics remain poorly understood. Here, we use a dynamical-systems approach to investigate how environmental conditions shape bloom dynamics from long-term time series of algal abundance and abiotic variables. Using the global modelling framework GPoM, we reconstruct nonlinear ordinary differential equation systems and identify candidate models capable of reproducing key features of the observed dynamics, including chaotic behaviour. Topological and statistical analyses reduce the candidate models to two systems that best capture the observed distributions and dynamical patterns. These models reveal distinct combinations of environmental conditions associated with bloom occurrence and amplitude, providing a framework for interpreting bloom trajectories and a potential basis for short-term prediction.

We develop a general geometrical framework for growth processes in which the growth rate depends exponentially on the number of active growth-supporting components. The framework separates the fraction of active components from the size of the accessible component space, allowing qualitatively different growth laws to arise from the same underlying construction. Bacterial population growth provides an example of a fixed-size component space controlled primarily by activation and deactivation, whereas examples such as human population growth or speciation provide an expanding component space in which interactions between individuals promote further growth. This motivates a Replicative geometrical law, dX/dt= X exp({beta} X^{gamma}), whose exponent characterizes the scaling of the effective growth-supporting component space and is tested against historical human-population regimes and eukaryotic lineage-through-time data. For innovation indicators such as patents and arXiv manuscript submissions, where an existing item need not itself reproduce, the same construction motivates a Direct Geometrical production law without the prefactor X. The framework therefore provides a common description of fixed and expanding component spaces while distinguishing Replicative from Direct Geometrical growth.

Known for their rapid water surface locomotion, whirligig beetles (Gyrinidae) are also fast divers. They can transition from the surface to underwater swimming in under 100 ms. Diving requires them to generate sufficient hydrodynamic torque to break the water surface and reorient the body. However, the mechanism by which whirligigs create the necessary pitch-down torque is unknown. Here, we show that their diving is achieved by bending the abdomen. The inertial hydrodynamic pressure on the flexed abdomen generates the torque necessary for surface breaking and body rotation. The torque is maximized at an abdominal angle of 45{degrees}, whereas for effective diving, whirligigs display a preferential abdominal angle of around 28{degrees}, closely aligning with the angle that maximizes the torque-to-drag ratio. This abdominal angle is optimized by balancing torque generation with the drag penalty. Our findings inspire a mechanism for surface-to-underwater transition, offering design principles for multimodal robots capable of locomotion across the air-water interface.

Spiders arrest their fall by abseiling on dragline silk, varying silk release rate to avoid excessive strains that could break the thread. To understand how this is achieved, we develop a minimal dynamical model of abseiling and formulate a time-optimal control problem subject to physical limits on silk strain and dragline spinning rate. We find that the optimal trajectory includes an interval over which the strain is constant, resulting in uniform deceleration consistent with published experiments. To test the generality of this strategy, we then turn to a distinct extreme motor behavior seen in jumping spiders: using a dragline anchored at the takeoff site to correct their orientation while jumping at speeds exceeding 100 body lengths per second. We show that the same constant-strain strategy suffices to counteract the initial body rotation and ensure an upright landing without significant loss of range. Together, our theoretical framework suggests how strain regulation may serve as a strategy for controlling extreme movements such as abseiling seen in many spiders and may be specialized to tasks such as aerial braking and reorientation in jumping spiders.

Historical plague mortality records provide unusually detailed temporal records of epidemic dynamics. However, the ecological processes underlying those dynamics are rarely observed. This creates a fundamental challenge for inference into ecological mechanisms, as multiple combinations of unobserved ecological parameters and initial conditions may produce similar mortality trajectories. We investigated this problem by examining the geometry of parameter compensation for models of the Barcelona plague outbreak of 1490. We fit a mechanistic rat-flea-human transmission model to daily human mortality data under three values of flea search efficiency, then characterized the resulting joint posterior distributions. We then used a phenomenological model in which human infection was driven by a time-varying external force of infection, allowing us to separate parameter dependencies arising from deterministic structure. The mechanistic model reproduced the broad temporal structure of the observed epidemic, but exhibited strong linear and nonlinear dependencies among transmission parameters and initial rat and flea conditions. Similar dependencies persisted in the phenomenological model, despite its substantially reduced complexity. However, fixing the timing of the external hazard substantially reduced these dependencies. These results demonstrate that detailed mortality records can constrain the temporal pattern of an epidemic without uniquely identifying the mechanisms that generated it. Rather than simply adding additional mortality records, improving mechanistic inference from historical epidemics may require recovering independent observations of upstream ecological processes, such as timing of rodent or flea dynamics as well as environmental drivers. Mechanistic models can therefore serve not only to reconstruct historical epidemics but also to identify which missing observations would most improve inference.

Buruli ulcer (BU) is a neglected tropical skin disease caused by infection with Mycobacterium ulcerans. The mode of its transmission into the skin has not been fully clarified, but it is believed to involve mechanical transfer, perhaps by insect bites. M. ulcerans secretes a toxin called mycolactone which causes necrosis of the skin and formation of large skin lesions. Understanding the development of these skin lesions over time is crucial to understand the pathogenesis of BU. We formulate a mathematical model to describe the space-time dynamics of M. ulcerans bacteria, mycolactone and skin cell densities as well as fibrin deposition. We investigate the relative contributions to ulcer formation of direct cytotoxicity of mycolactone and ischaemia following the loss of coagulation control. Our results are consistent with existing biological evidence that higher mycolactone concentrations accelerate cell death. Furthermore, the model indicates that the rate of wound expansion is primarily controlled by the diffusion and net birth rate of bacteria, which sets the overall scaling for the speed of wound expansion, along with two dimensionless constants that define the shape of the trajectory of wound edge position against time, namely the ratio of the mycolactone diffusion coefficient to the bacterial diffusion coefficient, and the ratio of the mycolactone degradation rate to the bacterial net birth rate. When mycolactone diffuses significantly faster than bacteria spread or when mycolactone degrades more slowly than bacteria reproduce, the initial development of the wound is faster. In contrast, the long-term wound expansion speed appears to be relatively unchanged across a wide range of the dimensionless parameters.

Fitting mathematical models of population dynamics to microbial time-series data allows us to estimate the ecological processes and interactions taking place in the microbiome. Repeated experiments of microbial systems yield replicates which slightly differ from each other. Some of this variability arises due to the fact that births and deaths occur at random. Most prior work focuses on fitting a deterministic mathematical model to the average across replicates. We use a stochastic model to fit the variability to the observed variability across replicates. Using a simulation-driven approach, we study the conditions under which our approach allows us to infer a larger fraction of ecological parameters correctly. We observe a substantial improvement in parameter inference. Lastly, our Bayesian approach not only allows us to incorporate prior information about the system, but also provides a distribution of parameters which conveys some idea of the uncertainty of the estimates.

Host immunity acquired through prior infections is a key factor shaping the spread of emerging variants. However, accounting for complete infection histories in mathematical models is challenging because the number of possible histories grows exponentially with the number of variants. Here, we present a flexible framework for modeling the co-circulation of multiple variants emerging at arbitrary times while incorporating cross-immunity and variant-specific characteristics. We introduce three tractable approximations: truncation, which retains only infection histories up to a given number of strains; one variant at a time, which is exact when variants circulation does not overlap; and two variants at a time, which is exact when at most two strains co-circulate. We systematically evaluate these approximations across a range of cross-immunity levels and inter-variant emergence delays. Truncation approximations perform best under short emergence delays and strong cross-immunity, whereas the other approaches outperform them at longer delays. When applied to SARS-CoV-2 variant dynamics, the second-order truncation and the two variants at a time approximation accurately capture incidence patterns, whereas the other approaches perform well only for specific subsets of variants. These results identify the parameter regimes in which simplifying assumptions remain valid, provide guidance for model selection, and highlight the key role of infection history in shaping epidemic trajectories.

A major goal in bottom-up synthetic biology is the construction of multicellular synthetic systems capable of coordinated and robust collective behaviours. However, robustness is often limited by noise and variability arising from increased molecular complexity. Whilst communication has been implemented in synthetic multi-cellular systems, the ability for communication to suppress cell free gene expression variability in populations of synthetic cells remain unexplored. To address this, we encapsulated the Lux and Las quorum sensing gene circuits in lipid vesicles under cell-free conditions to test the effect of communication on reducing cell-free gene expression variability across the population. Our results show that communication, limiting expression resources, and membrane surface effects can reduce gene expression variability. Resource limited Gillespie simulations for transcription and translation show that communication-mediated coupling reduces population-level expression noise under constrained and excess resource conditions. Together, our work provides simple strategies to reduce gene expression variability and thereby improve robustness in synthetic multicellular systems, an important criteria for the future applications of synthetic cells.

The emergence of treatment resistant mutants is one of the primary contributors to therapy failure in cancer. Evolution-based adaptive therapies aim to exploit the competition between sensitive and resistant cells to prolong therapeutic control. However, in the densely packed tissues of solid tumors, competition is not only determined by cell-intrinsic factors but also a result of local physical interactions between neighboring cells. The successful development of physics-informed adaptive therapies therefore hinges upon the bottom-up integration of tissue mechanics and experimentally grounded models of tumor evolution. Here, we show how population-level treatment response emerges from mechanical interactions between individual cells. Combining genetically tailored therapy-mimicry experiments in spatially structured 2D melanoma colonies with mechanistic, particle-based tissue simulations as digital twins in a real-to-sim-to-real framework, we reveal how treatment-induced changes in mechanical crowding shape a cascade of resistance expansion, spatial competitive escape, and consequent therapy failure. We then link these underlying physical mechanisms to population-level observables in an abstracted mathematical model for reinforcement-learning-based therapy optimization and a distilled treatment-decision-heuristic, which we test experimentally. Together, our results identify mechanical crowding and escape as governing principles of resistance-mediated therapy failure in solid tumours and provide a basis for mechanistically informed, evolution-based treatment strategies.

Species at the highest trophic levels are generally the least resilient to perturbations, yet the factors that determine their persistence remain poorly understood. We ask how the intrinsic self-regulation of a basal host governs the persistence of the trophic levels above it. To address this question, we use a minimal discrete-time model of a host-parasitoid-hyperparasitoid system motivated by the biology of trigonalid wasps, obligate hyperparasitoids that complete development only when their caterpillar host is also attacked by a primary parasitoid. We compare two forms of host self-regulation in the absence of natural enemies: scramble (Ricker) fecundity, appropriate when the parasitoid is the host's principal enemy and which recovers the classical model of Beddington and Hammond as a special case; and contest (Beverton-Holt) fecundity, appropriate when the host is held in check by many additional enemies. Using numerical continuation, Lyapunov exponents, and direct computation of the model's trajectories, we find that the form of host regulation is decisive for persistence at the top trophic level. Under scramble host dynamics, coexistence is extremely fragile with respect to both parameter values and initial conditions. Under contest host dynamics, coexistence is robust across broad parameter ranges and loses stability only gradually, through a supercritical Neimark--Sacker bifurcation and a subsequent torus-doubling route to chaos. Therefore, in at least some trophic networks, the resilience of the highest levels can be determined at the base of the network.